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FROM TRANSIENT-RESPONSE DATA 


By Marvin Shinbrot 


SUMMARY 


A general theory of the so-called "equations -of -motion” methods for 
the analysis of linear dynamical systems is developed first. It is then 
shown that when viewed from this general point of vantage, all of these 
linear methods can be extended in a straightforward manner to apply to 
the analysis of nonlinear systems. In addition, through use of this 
theory, a new method is derived. It is essentially a variation of the 
well-known "Fourier transform" method for the analysis of linear systems 
but possesses certain advantages over previous methods. Application and 
effectiveness of this method are demonstrated by three examples, two of 
which are nonlinear - one highly so - and the third being of the fourth 
order . 


INTRODUCTION 


It has often been suggested (e.g., in ref. l) that nonlinearities 
which are ignored in the classical theory of the equations of motion of 
an aircraft may be responsible for certain unusual phenomena which have 
been observed in flights of modern high-speed airplanes and missiles. 
Consequently, it seems desirable to develop methods for the analysis of 
such nonlinear systems - methods which allow the calculation from measured 
transient-response data of the nonlinear stability characteristics as well 
as the classical linear stability derivatives of the aircraft. Several 
such methods are described in reference 2, the principal one consisting 
of a generalization of the so-called "derivative method" which was orig- 
inally devised for use with linear systems (cf. ref. 3)* However, the 
methods described in reference 2 leave something to be desired from both 
points of view of accuracy of the results and lengthiness of the calcu- 
lations. In addition, application of these methods requires, in all but 
the simplest cases, the previous evaluation by some means of those sta- 
bility characteristics which are linear. In view of these shortcomings, 
an attempt has been made in the present study to find simpler, more accu- 
rate, and more general procedures. The problem is attacked by first exam- 
ining several well-known methods for the analysis of simple linear sys terns 
and then modifying them as necessary to allow their application to more 
general systems. 
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Many methods for the analysis of linear systems have been proposed 
in the past (see, e.g 0 , refs. 3 and 4). In reference 5, these methods 
have been classified under two main heads: "equations -of -motion” methods 

and "response-curve-fitting" methods, the former title including the 
derivative method and what have been called the Laplace transform and the 
Fourier transform methods (r*ef. 3 ), and the latter consisting of such 
methods as Prony f s (refs. 3 and 6) and the techniques of reference 4. 
Since the response -curve -fitting methods involve the explicit— solution 
of the equations of motion in terms of the physical parameters of the 
system at hand, they do not seem suitable for use with nonlinear systems. 
Hence, we shall-be concerned solely with the equations -of -motion methods. 

Each of these methods has been considered in the literature as an 
independent, entity; apparently, no attempt has ever been made to subsume 
all of them under a single general theory. For the purposes of the pre- 
sent study, such a theory would be desirable since it seems reasonable to 
expect-first that when viewed from a more general point of view, a gen- 
eralization of the methods to nonlinear systems might appear; and second 
that once such a theory is known, it might be possible to develop new 
methods, superior in certain respects to the old ones from which the 
theory sprang. 

In accordance with this plan, the paper begins with a short - presen- 
tation of the three best known of the equations -of -motion methods. These 
methods are examined from a new point of view which is then shown to lead 
to the general theory for linear systems. The further extension to non- 
linear systems is considered next. Based on the general theory, develop- 
ment of a new method for data analysis follows. Finally, some examples 
of the application of this recommended method are given. 


SYMBOLS 


a 

a 

b 


Co 


c i 

& 

% 


angle of attack, radians 

da 

dt 

Ja Mj | 

mV ‘ I y ' Iy 

% I« Mq 

Iy + mVIy 

L& 

mV 

elevator deflection, radians 
pitching moment of inertia, slug-ft 2 
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k 

L 

% 

1$ 


- ife _ Ja !k 

“ Iy “ mV Iy 

linear lift force, lb 


8l 

8a 

8l 

85 


L(a,) 

m 

M 

Ma 

Ma 

Mg 

«<i 


nonlinear lift force, lb 
mass of aircraft, slugs 
linear pitching moment, ft-lb 

8m 

8 a 

8m 

8 a 

8m 

85 

8m 

8q 


M(a) nonlinear pitching moment, ft-lb 

q pitching velocity, radians/sec 

t time, sec 

V velocity of aircraft, ft/sec 

Other symbols will be defined as they are introduced. 


ANALYSIS 

Some Equations -of -Mot ion Methods for Linear Systems 


In the presentation of the general theory of equations -of -mot ion 
methods, it is desirable to have some examples of these methods set down 
for examination. The three best known of these methods are briefly 
described below; for a more detailed discussion of them, see reference 3 . 
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As a concrete example, let us consider an airplane operating under 
conditions where the stability characteristics are 'effectively linear, 
so that, as in reference 3; the equations of its longitudinal motion can 
be written r 


- ccl^ - amV + qmV = 5L§ 

- ogMq^ - aWa - qMq + qly = SMg 

It will be convenient to eliminate q and so to write these equations as 
a(t) + ba(t) + kn(t) = C 0 5(t) + C^Ctr) (l) 


where dots denote differentiation with respect to time t. It is assumed 
that time histories of a(t) and 6(t) are available from flight records 
and that one wishes to calculate the constants b, k, C 0 , and C 1 . For 
simplicity, it will further be assumed that the constant b is positive 
and that &(t) is a pulse-, so that 8(t) is zero for t sufficiently 
large. These restrictions will be removed farther on. 


The , derivative method . - In order to apply the -derivative method, It 
is necessary first to differentiate the given records to obtain the 
derivatives a(t), a(t), and S(t). Then, for fixed ty equation (l) may 
be considered as an equation for the constants b, k, C G , and C x . By 
letting t vary, a set of such equations can be obtained, which can in 
turn be solved by least squares (see reference 6, page 210, where, however, 
the integral rather than the sum of the squares of the errors has been 
minimized) for the desired constants. By this procedure, the following 
equations are obtained: 


r*°°. p 00 p°° * p m ^ r>°° 

b J a^dt + kj oodt - C 0 o5dt - C ± J a Bdt-=— = J a, adt 


poo ' pCO pOO pCO * pCO 

b / cxadt + k / a 2 dt- - C 0 / abdt - Cj/ / ocSdt = - / aadt 

'-'o ^0 

p°° poo pOO pOO t pOO 

- b / Sadt - k / Badt + C Q / S 2 dt + C 1 / 88dt = / Sadt 

Ufi '-'n '-'n 


>( 2 ) 


poo ^ ^ poo ^ pec _ poo . poo ^ 

- b j 8 adt -k J 8adt + C Q / 86dt- + C 1 j S 2 dt = / 8 adt J 

Equations (2) can be solved for the desired parameters b, k, C Q , and C^. 
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The Laplace transform method .- Letting A(p) and A(p) denote the 
Laplace transforms of <x(t) and S(t), respectively, bo that 


A(p) = P e" p ^ a(t)dt 
J o 

A(p) = e _ P t 6 (t)dt 
J o 

it follows that if a,(t) and S(t) are related by equation (l), then 


(3) 


(p 2 + bp + k) A(p) = (C x p + C 0 ) A(p) 


(4) 


(ref. 7 ). In writing down equation (4), it has been assumed for simplic- 
ity that a(o) = a( 0 ) = 5(0) = 0 ; this restriction is inessential and will 
be removed later on. For any value of p, equation (4) is an equation 
in b, k, Co, and C 2 . After finding the Laplace transforms of a(t) and 
5(t) for several such values of p, say for p = p 1 , p 2 , . . ., pjjj, the 
corresponding equations (4) can be set up and solved by least squares to 
obtain 


> (5) 


h X Pi 2A - 2 (Pi) + k X p i A2 ( p i^ ' C ° X PlA ^ Pl ^ ‘ 

^ Pi 2 A(p i ) A(Pi) = ' X p i SA2 ( p i^ 
b PiA 2 ( Pi ) + k X a2 (p±) - Co X A (Pi) A( Pl ) - 

c i X p i A ^ p l^ = " X p i 2A2 ( p i) 

- b ^ Pi A (Pl) A(Pi) - k X A (Pi) A (Pl) + C 0 A 2 (pi) + 

C l X P l Aa ^ P i^ = X Pl - 2A ^ Pi ^ A ^ Pl ^ 

- t> ^ Pi 2 A( Pl ) A( Pi ) - k X Pi A (P±) A (Pl) + c o ^ PiA a (Pi) ■ 

C x ^ Pi^CPi) = X p i 3A ^ p i^ A (Pi) 

where all sums are over the range 1, . . . , N of the index i. These 
equations then can be solved for the desired constants. 


The Fourier transform method .- The Fourier transform method proceeds 
in much the same way as did the Laplace transform method. 
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Defining 


A(ico) - P a(t) dt 
J 0 

A(iw) = f~e- iat 8(t) dt 
(i 2 = - l), it follows that- if a(0) = <x(0) = &(0), then 


( 6 ) 


[(ico) 2 + icob + kj A(ico) = [C Q + iwC-J A(i<o) ( 7 ) 


ThiB equation can be written as two real equations by setting 


A(ico) = C(u>) - iS(to) 

A(ico) = r(u) - i 2 (to) 

Putting equations (6) and (8) together, we see this__means that 

/ CO poo 

a(t) cos cotdt S(co) = / ac(t) sin cotdt - 


T(w) = P 6(t) cos 
Jo 

Breaking equation ( 7 ) into its real and imaginary parts gives 
co S(co)b + c(to)k - r(a))C 0 - co Z(co) C 1 = c^CCco) 

to C(co)b - S(co)k + E (to)C 0 - co r(to)c ;L = - w 2 S(co) 


(Otdt 


z(co) = J' 6(t) sin cotdt 


( 8 ) 


After C(co) , S(co), F(co), and Z(co) have been evaluated for several values 
of co, substitution into these equations yields a number of-^equations 
for b, k, C 0 , and C- L which can be solved by least squares. These equa- 
tions, corresponding to (2) and (5)> are exceedingly complicated and will 
not be reproduced here. 
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The Common Feature of These Methods and the First Generalization 


Each of these methods is usually derived as in the preceding section, 
through use of a certain specialized concept. Thus, the derivative method 
is based on the fact that for fixed t, equation (l) is linear in b, k, 
C G , and C-l, so that a set of simultaneous linear equations for these 
parameters can be -obtained by varying t. The Laplace and Fourier trans- 
form methods stem from the theory of these operators. Thus, if a(t) 
and S(t) are related by the linear differential equation (l), their 
Laplace transforms are related by equation ( 4 ). However, equation ( 4 ) 
is linear in the parameters, and so by writing this equation for several 
values of p, one can solve the resulting equations by least squares for 
the desired constants. These derivations tend to obscure the common 
idea which can be shown to lie behind all the methods. This difficulty 
can be overcome if these particular derivations are forgotten and if 
attention is fixed entirely on the formal processes whereby the final 
least squares equations are obtained. With this in mind, let us recon- 
sider the three methods. 

The derivative method . - Equations (2) for the derivative method are 
formally obtained by 

(1) Multiplying equation (l),by the four functions 

a(t), a(t), S(t), and S(t) one at a time and 

(2) Integrating the results from zero to infinity. 

We should forget for the moment the interpretation of this procedure as 
the solution of equation (l) by least squares, and simply keep the process 
of multiplication and integration in mind. 

The Laplace transform method .- A similar process can be described 
for the Laplace transform method. Choosing N( > 4 ) positive numbers pj_, 

(1) Equation (l) can be multiplied by the functions 

e-Pit, i = 1, . . ., N, and 

(2) The results can be integrated from zero to infinity. 

If the resulting equations are solved by least squares, precisely equa- 
tions ( 5 ) for the determination of the parameters by the Laplace trans- 
form method are obtained, provided that in step ( 2 ) any integrals which 
arise involving derivatives of a(t) and b(t) are integrated by parts 
to eliminate these derivatives. 

It should be noted that although there appears to be one more step 
here than there was in applying the derivative method - notably an addi- 
tional least squares following step (2) - this addition is more apparent 
than real, since it is necessary to apply a least squares process here 
merely because N (which is generally greater than four) equations are 
obtained for the four parameters, while in the derivative method exactly 
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four such equations were obtained. Thus, the least squares step could 
be eliminated by choosing N = 4 (of course , such a choice is not really 
practical), or, alternatively, such a step could be added to the deriv- 
ative method by solving the equations resulting from step (2) in that 
method by least squares instead of in the ordinary -way. Of course, such 
a step would only be made to make the two methods so far described 
formally more similar; in practice, it would not be performed. 

The Fourier transform method .- Finally, 

(1) If equation (l) is multiplied by cos tot and sin cot 

for several values of co, and 

(2) If the results are integrated from zero to infinity 

(as in the Laplace transform method, integrating by parts 
to eliminate explicit dependence on the derivatives of- a, 
and 8), 

one obtains a set of equations identical with those obtained from the 
Fourier transform method. •— 

The general method for linear systems .- The general development of~ 
equations -of -mot ion methods is now manifest. One takes the equations of 
motion for the physical system under consideration - for definiteness, 
say equation (l) - and 

(1) Multiplies them by N arbitrary (but sufficiently 

smooth) functions y v (t) . 

(2) The resulting equations are then integrated between 

two definite limits, say, zero and T, 

In the three methods described above, T = co, but this is not essential. 

In order to avoid some complications initially, we shall continue to 
integrate over this infinite interval; this restriction will subsequently' 
be removed, however, and T will be allowed to have finite values. In 
the case of equation (l), the process just described leads to N equa- 
tions of the form 


/ CO ^ poo pCO 

y v (t)a(t)dt + k / y v (t)a(t)dt - C 0 / y v (t)5(t)dt-=- 
U J o J o 

/ CO 4 pOO 

y v (t)8(t)dt-= - J y v (t-)a(t-)dt^— V = 1, . . . , N (9) 

o o 


It is possible that the functions y v (t) depend on a(t) or S(t) as, for 
example, in the derivative method; in such cases, equations ( 9 ) can be 
considered as N equations which are to be solved by least squares for 
the desired parameters. Of course, this process requires the calculation 
of the derivatives a(t), a(t), and 6(t). On the other hand, if the 
functions y v (t) are explicitly independent of a and 5, as Is the case 
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in the Laplace and Fourier transform methods', the following formulas, 
obtained by integrating by parts, are used: 

P y v (t)a(t)dt = - y v (0)a(o) - P y v (t)a(t)dt 

J o J o 

/ oo pco 

y v (t)a(t)dt = - y v (0)a(0) + y v (0)a(o) + J y v (t)a(t)dt 

/ OO ^ ptO 

y- v (t)8(t)dt = - y v (0)8(0) - / y y (t)8(t)dt 

u J o 

Substitution into equation ( 9 ) gives 

- b j” yv(o) a (0) + P y v (t)a(t)dtl + k P y v (t)a(t)dt - 

C 0 >P y v (t)B(t)dt + C-lT y v (0)S(0) + P y v (t)&(t)dtl = y v (0)d(0) - 

o L <-/ 0 -* 

/ CO 

y v (t)a(t)dt, V = 1, 2, . . N (10) 


Equations (10) are N equations in b, k, C 0 , and C x . If N > 4, they 
may be solved by least squares for these parameters. 

The choice of the functions y^(t) defines which equations -of -motion 
method is being used. For this reason, these functions will be referred 
to as the "method functions . ,T 


Generalization to Nonlinear Systems 


In this section, all considerations -will refer to the equations of 
longitudinal motion of an aircraft. It will be seen that the method 
actually is applicable to a far wider class of equations - in particular 
to the equations of lateral motion of an aircraft, including, if it is 
so desired, cross -coupling terms. The special analysis presented here 
can be extended to other problems, the only real restriction being that 
for practical reasons too many parameters cannot be handled at once. 

The following equations, which Involve assumptions of constant air- 
speed, smallness of certain quantities, etc., are often used to describe 
the motions of an aircraft which has linear stability characteristics: 

- al^ - amV + qmV = 5L$ 

- a^a - - qMg + qly = SMg 
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If, on the other hand. It is assumed that the lift and moment functions 
are nonlinear in a, these equations become' 


- L(cg) - amV + qmV = SLg 
- M(a) - dl& - qMq + ql y = 61% 


We shall assume that over the range of interest, approximations of the 
following form are valid: 


L(a) = L-jCt, + I^a 2 + . . . + L -^ cl 1 
M(a) = M x a + Mga 2 + . . . + M n a n 


( 12 ) 


where the "coefficients L^, M-^ are constant. Only the first three terms 
of this series will be retained in the present- analysis since this three- 
term approximation usually balances very well the opposing requirements 
of simplicity and adequate representation of the aerodynamic parameters. 
If more terms are found to be necessary in a particular problem they can, 
of course, be added. It should be noted that this three-term approxima- 
tion is still fairly general, even retaining the possibility of asymmetry 
in the nonlinearities . 

By use of the approximations (12) with l - n = 3> equations (ll) 
can be written 


(13) 


The generalization of the methods of the preceding section to such non- 
linear systems proceeds in the obvious way. First, multiply each of 
equations ( 13 ) by N method functions y v (t) which have been selected 
as suitable. This operation is followed by integration of the resulting 
2N equations. If records of-the derivatives q and a are not available, 
one then eliminates the terms involving these derivatives by integration 
by parts. The result is N equations In L x> L 2 , L 3 , and Lg and N equa- 
tions in M x , M 2 , Mg, M&, and Mg. If N > these two sets of equations 
can' be solved by least squares for the parameters. 

Of course, the pitching velocity q can be eliminated from equa- 
tions (ll) to yield the single equation t 


l i L 3 . 1$ 

a — or - — a°-a+q = — 5 

mV mV mV * mV 


% . Mq . Mg 

— a - q + q = 5 

x y x y x y 
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„ [m,M a L'(co) 1 M(a) Mq L(a) Lg * 

■ l Iy Iy " mV J Iy ’ Iy mV V*y “V ° ~ mV ° 

where 

L*(a) - £ L(a) 

With the approximations (12), this leads to the following generalization 
of equation (l): 


a + (b Q + + b 2 a, 2 )a + (ko + k-jCt, + k^a^a = C 0 8 + 


where 



(14) 


(15) 


By applying the method described to equation (1*4-), the constants b 0 j b^ 
b 2 , ko, k-L, k 2 can be calculated. 


Solution for the parameters in equations of the general form of 
equation (l^) is of interest to workers in many fields. In addition, 
the evaluation, from these parameters, of the stability constants occur- 
ring on the right sides of equations ( 15 ) is of considerable interest to 
aeronautical engineers, and so this problem will be considered in further 
detail. One cannot, in general, isolate the constants and Mi to 
obtain, from values of the bi and the kj_, the nonlinear functions L(a) 
and M(a) . For this reason, it is ordinarily best to apply the method 
directly to equations ( 13 ) rather than to equation (l4), provided that 
records of both q(t) and a(t) are available. 

On the other hand, there is one case of interest when equation (l^) 
can be used directly and measurements of q(t) are not needed. It some- 
times happens that while the pitching moment M(a) is nonlinear, the lift 
L(a) can still be successfully approximated by a linear function. In 
this case, we have the approximations 

L(cc) = IfcOt 

M(a) = M x a + M^x 2 + M©a? 


( 16 ) 
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and so equation (l 4 =) becomes 


a + b& + (k D + k 2 a + k ga 2 )a = C G & + C 18 


(17) 


where 


k o 

*1 


Mq IaV 

\ ly 


(18) 


k - Ms 

ly 

Thus in this case, M(a) can, except for the term & ~ occurring in the 

Iy mV 

expression for ko, be obtained from an analysis of- equation (17)- At 
high speeds, however, this term is small and its effect on the curve o£_ 
M(a) versus a, can be neglected. Thus, the expression for ISq in 
formulas (l 8 ) can be replaced by the expression 


" - $ U9) 

and in this case, the nonlinear moment M(a) is completely determined by 
the knowledge of ko, k a , kg, and, of course, I y . 


Choice of the Method Functions 


Up to this point in the. general discussion, the method functions 
y v (t) have been to a very great extent arbitrary, having to satisfy only 
certain weak smoothness conditions. In this section, the possiblity of 
developing new and perhaps improved methods for data analysis by means 
of a particular choice of the method functions will be explored. 

Previous experience, consisting in part of unpublished analyses per- 
formed at Ames Aeronautical Laboratory, have indicated that the Fourier 
transform method generally results in greater accuracy than either the 
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derivative or the Laplace transform methods. 1 For this reason, the 
method to be discussed -will be a variation of the Fourier transform 
method. 

Let us consider first some evident shortcomings in the Fourier 
transform method. These defects will offer definite goals to be held 
in mind in the development of a new method. First, there is a weakness 
in the Fourier transform method in that all integrations proceed over 
the interval from zero to infinity (cf . eq. (10)). This causes diffi- 
culties in any example in which cc(t) and 8(t) do not approach zero so 
rapidly that their integrals exist. Thus, referring, for example, to 
equation (l), if b is not positive or if 5(t) does not approach zero 
quickly enough, the method cannot be applied straightforwardly . Further- 
more, even if b >0 and 5(t) — > 0 in such a way that /°°a(t)dt exists, 
the experimental record often is not long enough for this integral to be 
accurately calculable when the system is so lightly damped (b small) 
that sizable oscillations persist even to the end of the run. One device 
which is sometimes used to overcome this difficulty is equivalent to a 
change in the method functions. Instead of the functions sin tot and 
cos tot, the functions e - ^ sin tot and cos tot, with some fixed 

constant are used. However, this leads to the same objection that 
was voiced in footnote 1 for the derivative and Laplace transform 
methods, notably that the method functions approach zero. Other tricks 
for dealing with such deficiencies iij. the Fourier transform method can 
be evolved, but, rather than develop new devices for each special case, 
it appears wiser to construct a generally applicable method in which 
these difficulties never arise - that is, one in which the integration 
proceeds only over a finite interval. 

The second defect which we shall consider becomes clear from an 
inspection of equation (10), with the functions yy(t) of the form 
sin tot and cos tot for certain values of to. Referring to equation (10), 
it is easily seen that one point, namely the point t = 0, is weighted 
very heavily because of the occurrence of the quantities a(0), a(0), 
and 6(0). It should be noted that not only are the values of a,(t) and 
8(t) at one point relied on to this great extent, but that even the 
relatively inaccurate value of the derivative of a(t) at that point is 
weighted. Thus, advantages in accuracy might be expected to accrue if 
these terms were eliminated. 

Since the method for overcoming this second deficiency in the Fourier 
transform method will also be used in treating the first, the problem of 
eliminating dependence on the initial values will be discussed now. To 
this end, consider equation (10). We begin with the Fourier transform 

1 It would seem that no rational explanation for this conclusion has 
heretofore been offered. However, the theory described herein appears to 
afford such an explanation. A long and rather tedious analysis based on 
this theory has indicated that the failure of both the derivative and the 
Laplace transform methods is due in large part to the fact that the asso- 
ciated method functions approach zero very rapidly as time progresses . 
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method , that is, with method functions yv(t) of the forms sin wt 
and cos wt. Noting that mosl^of the terms which depend on the initial 
values at 1 r~^ 0 are multiplied by y v (0), it is seen that a choice of 
method functions such that y v (o) is zero for all V - 1, 2, . . . , N 
represents a step in the right direction. Such a choice is easy to 
make, since it is only necessary to eliminate those method functions 
which have the form cos wt; then, we may write 


y^(t) = sin w v t, V = 1, 2> . . . , N 


(20) 


This does not entirely eliminate the dependence on the initial conditions, 
however, as the term y v (0)a(0) remains in equation (10). If this term 
can also be removed, the second weakness in the Fourier transform method 
will have been entirely corrected. This will clearly be the case if 
y v (0) as well as y v (o) is zero for all V = 1, . . ., N. A possible 
choice of the method functions for which this is so, a choice which 
still retains the advantages of the favored Fourier transform method, is 
the following: 


y v (t) = sirf^Wyt = 


1 - cos 20 )yt 


V = 1, 


N 


(21a) 


With this choice of the method functions, equation (lO) becomes 

/ » ptO ptO 

a '(t)y v (t-)dt + k / a(t)y v (t)dt - C Q / B(t)y v (t)dt 
u a o a o 

/ oo B p» 

S(t)y v (t)dt = - /• a(t)y v (t)dt 


( 22 ) 


The method functions (21a) would be used for systems satisfying dif- 
ferential equations, like (l), which involve derivatives of the second 
order. More generally, and for the same reasons, if the highest order 
derivative occurring in_any of the equations of motion of a system is 
the nth, the following method functions are suggested: 


y v (t) = sin 11 w v t, V =s 1, . . . , N 


(21b) 


Thus, in the case of the two -degrees -of -freedom system described by 
equations ( 13 ), there is no point in using formula (21a); the simpler 
method functions (20) may as well be used. 

As for the first of the weaknesses in the Fourier transform method, 
that which is due to integration over an infinite interval, it would 
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appear at first glance to be easily disposed of by merely choosing some 
finite, positive number T and integrating over the interval from zero 
to T. This, however, introduces another difficulty. Suppose yyi(t) is 
given by equation (21a), so that y v (0) = yv(°) = Returning to the 
derivation of equation (10), it may be seen that if one integrates only 
over the interval 0<t<T, all the good which has been achieved by 
eliminating dependence on the point t = 0 is obviated by certain terms 
which arise - terms of the form y v (T)a(T), y v (T)a(T), and y v (T)&(T). 
Thus, the heavy dependence on the initial conditions is replaced by a 
dependence on the final conditions. Of course, the same approach as 
was used for eliminating the initial conditions can be used again - that 
is, the method functions can be chosen in such a way that 
y v (T) = y V ( T ) = Gne Possibility, naturally, is to choose the frequen- 
cies o)y such that 


sin ojyT = 0, 


V = 1, 


N 


Thus , we can set 


u v = V = 1, . . ., N 


This choice of the frequencies (corresponding to the method functions 
given by equation (21b)) leads to an elegant method which gives satis- 
factory results in certain cases. On the other hand, the difference 

^ between two successive frequencies is too large to define the 
requency response" (to use loosely the terminology of the Fourier 
transform) of some examples adequately. For this reason, we should like 
to be able to choose the frequencies as follows: 



V = 1, 


N 


(23) 


Of course, this means that for frequencies having an odd subscript, y-y(T) 
will be different from zero. To overcome this difficulty, the following 
choice of the method functions can be made: If the highest derivative 

occurring in the equations of motion is the nth, define the method func- 
tions by the formulas 




= sin 11 u 2M t 


f sin 11 t, 0 < t < 




l 0 , 


_2 M_ 

2p+l 


2u 

2p+l 


T < t < T 


(2fc) 
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where the frequencies Ofy are given by equation (23) and T is the 
length of the run. What has been done in choosing these method functions 
is the following: Those method functions which are such that yv(T) is 
zero (i.e., those method functions which have even subscripts), have been 
left unaltered in the form given by equation ( 21 b). The remaining method 
functions have had the last quarter cycle, during which they would nor- 
mally have varied from 0 to ±1, chopped off, so that they are identically 
zero over part— of the interval 0 < t < T. The claim is not made that- 
this is the best of all possible choices for the method functions. There 
is certainly no reason for such an assertion, particularly in view of the 
fact that a certain amount of the data is not used by each of the odd- 
numbered method functions. However, this amount is small after all and 
no datum is completely discarded, since the even-numbered method functions 
use all the data. That these method functions do seem to be adequate is 
indicated by the results obtained in the examples given below. 

Before- proceeding to the examples, one further change, imposed in 
the interest of simplicity in the computations, will be made in the method 
functions. For the odd-numbered method functions, certain complications 
arise in the computations due to the fact that the point (2u/2u+l)T at 
which the function is cut-off may not— coincide— with a point at which the 
data are tabulated. For this reason, values of even-numbered frequencies 
will be chosen in accordance with equation ( 23 ). The odd-numbered fre- 
quencies, on. the other hand, will be changed (by as little as possible, 
to be sure) in such a way that the following condition is satisfied. Let 1 
the data be tabulated at the points t = t Q , t^, . . ., tg, where 
^k+i ~ tk = constant. The odd-numbered frequencies are then assumed to 
be chosen such that' 


sin to^vH-itqr, _ 0 


where t a is that one of the two tabulation instants closest to the time 
( 2 u/ 2 u+l)T having an even subscript. The -reason for this last condition 
is merely that it is convenient for a numerical integration procedure 
(such as the one in the Appendix) to have an even number of intervals 
(tk, tk+i); thus, Simpson* s rule, for example, calls for an even number 
of intervals in its application. ~~ 

For many of the experiments performed on airplanes and missiles, the 
run is 2 seconds long and the time between tabulated points is 0.05 sec- 
ond, so that the tabulation times t a are at 0, 0 . 755 * 0.10, . . ., 2.00 
seconds. For these values of T and At, the rule given above for the 
frequencies cjy is climaxed by the following table: 


V 

2 

3 

1* 

5 

6 

7 

8 

9 

10 

11 

12 

13 

Ik 

15 

16 


* 

2 

10* 

13 

it 

5* 

T- 

3* 

2 

30* 

17 

2* 

20* 

9 

5 * 

2 

25* 

9 

3* 

10* 

3 

7* 

2 

70* 

19 

ij-Tt 
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The frequency has heen omitted from the table since, according to 

the rule for the determination of the method functions, y-j/t) = 0 
(see eq. (24), ff.) and offers no information. A 2-second interval will 
be taken as standard in this report, and all computations will be based 
on such an interval. If one has a data run T seconds in length, it 
can of course be brought into this standard form by a preliminary sub- 
stitution of the form 



and we always assume such a transformation has been made. 

Thus, finally, the method may be summarized as follows: Select N 

frequencies in accordance with the rule given in the preceding paragraph. 
(The number N is chosen, as in the Fourier transform method, large 
enough to cover the frequency range of interest in the particular problem 
usually, N = 16 is adequate.) Multiply the equations of motion by the 
method functions (24) and integrate the resulting equations from zero 
to T, where T denotes the length of the data run. Eliminating all 
explicit dependence of these equations on derivatives of the ‘data by 
successive integrations by parts results in N linear simultaneous equa- 
tions for the parameters . The coefficients in these equations are all 
integrals involving the recorded data; after these have been evaluated 
by some means, 2 the equations can.be solved by least squares for the 
desired parameters . 


EXAMPLES 


Three example problems will be solved in order to demonstrate the 
effectiveness of the proposed analysis method and to illustrate asso- 
ciated computing techniques. An effort was made to select examples 
representative of problems which often occur In aircraft-response flight 
testing and which have not been handled adequately by other known analy- 
sis methods. These examples have been simplified in some respects, not 
because of fundamental limitations of the method, but in order to avoid 
obscuring the essentials of the method and of the related computing 
techniques. 

Although limited use was made of automatic digital computing machin- 
ery in the following analyses, It did not appear worthwhile to mechanize 
complete calculation procedures for these isolated illustrative examples . 
However, the method appears to be well suited to such mechanization. 


2 See, for example, the Appendix, where a technique well suited for 
the type of integrations needed for this method is described. 
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Example I 


The first example concerns the longitudinal response of a hypothet- 
ical missile for which it is assumed that the lift varies linearly with 
angle of attack, although the pitching moment does ;hot. The velocity of_ 
the missile is assumed to be sufficiently high for_the expression (l8) 
for k Q to be simplified to 


ko = - 


Ml 


so that equation (IT) can be used to determine the lift and moment char- 
acteristics from a transient-response. A pulse response of a Bystem 
described by equation ( 17 ) was- obtained from a Reeves Electronic Analogue 
Computer and. it was decided to consider this response as given data to be 
analyzed. The moment, of inertia of the. missile was chosen to be 100 slug- 
feet . The nonlinear moment curve M(a) versus a which was used to 
obtain the data is shown as the solid curve in figure 1. The linear sta- 
bility derivatives were chosen in such a way that the- damping parameter 
is given by 


b = 2 (25) 

In order to simplify the presentation, it was decided that free 
oscillations alone would be analyzed to determine only the constants 
occurring on the left-hand side of equation ( 17 ). Thus, for the data 
which will be analyzed, S(t) - 0, and equation ( 17 ) becomes 


a + ba, + (ko + k^po + kgOp)a = 0 ( 26 ) 

A plot of the a(t) "data” is given in figure 2, and this information is 
listed in table I. — 

Since equation (26) is of the second order, we shall choose, accord- 
ing to the rule given earlier, the method, functions (2*4-) with n = 2; 
since the run is 2 seconds long and the time interval between data points 
is 0.05 second, the frequencies are chosen as in the table on 

page 16 . 

Integrating factors r u (y v ) corresponding to each function yy(t) 
and its first two derivatives are tabulated in columns 6 through JO of 
table I. (As discussed in the Appendix, the T n (yv) are. numbers chosen 
such that for any integrable function x(t), the sum 

2^ x{t n )r n (y v ) 

n 
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is an approximation to the integral of x(t)y v (t).) Accordingly, the 
sums displayed beneath table I are the integrals needed for the reduc- 
tion of the data, divided by the factor At = 0.05. Since this factor 
occurs homogeneously throughout the equations from ■which the parameters 
are to be found (i.e., the generalization of eqs. ( 22 ) to the nonlinear 
eq. at hand), it can be divided out of these equations, and the sums can 
be used directly without first multiplying them by At. 

In accordance with the method as it has been described, the equa- 
tions which are to be solved by least squares have the form 


— — [- b P a(t)y v (t)dt + ko f cx(t)y v (t)dt + k x P a 2 (t)y v (t)dt 

°.°5 L J D J 0 J o ■ 

*2 f a 3 (t)y v (t)atj = - a(t)y v (t)dt 


The sums below table I are needed for the evaluation of the coefficients 
of b, ko, k-L, and kg in the above equation. These sums are again 
liBted in table II, and the coefficients in the last equation are set 
down as columns 4, 5.? 8, 7 3 and 9 of table II. The sums displayed beneath 
table II are needed for the final least squares step of the solution. 

Using these sums, it can be seen that the following equations are to be 
solved for the parameters: 


24.1669 b - 1.05270 ko - 0.336760 k x - 0.0140294 kg = - 7.01790 
- 1.05270 b + 0.885340 ko + 0.144127 k-L + 0.00906192 kg = 45.5226 
- 0.336760 b + 0.144127 ko + 0.0374138 k x + O.OOI 89656 kg = 7.00244 
- 0.0140294 b + 0.00906192 ko + O.OOI 89656 k x + 0.000108975 kg = 0 . 459829 

(27) 

It may seem odd to some that four significant figures' have been used in 
table II for the values of the integrals and six significant figures for 
the coefficients in equations ( 27 )> while the test data are not given to 
more than three significant figures. The reason for carrying more sig- 
nificant figures in the computations than there are in the data is to 
avoid eventual loss of accuracy due to round-off errors and other errors 
of a similar type. This procedure of carrying a few more (fictitious) 
figures than the data supply is usually necessary in order to retain even 
the basic information which is in the data. 
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Solving equations (27) gives 

b - 1.95 
ko = 50*4 
= - 30.5 
kg = 806 


( 28 ) 


The only parameter -whose numerical value is given and -which can "be 
immediately checked is the damping parameter. Comparing the values 
of b given in equations ( 25 ) and ( 28 ), we see that it- has been found 
with an error of 2.5 percent. The constants ko, k 1 , kg cannot~be 
checked directly; however , the calculated pitching -moment- curve 


M(cc) = - Iy(koa + k^ 2 + l^a, 3 ) 

= - 5040 a + 3050 a 2 - 80600 a 3 


can be plotted and compared with the true curve from which we started. 
This has been done in figure 1 from which it can be Been that the error 
at the least accurate point is less than 3 percent; 


It should be noted that the values of T n (y v ) given in table I can 
be used to solve any problem of the type considered here which depends 
on a second-order, differential equation or on a system of-such equations. 
If-J:he data run is 2 seconds long, Itr^Ls only necessary to insert the 
data in. table I in place of the data used In this example-and proceed 
as we have just done-. As mentioned earlier, if the data run is more or 
less than 2 seconds long, it is only necessary to make a preliminary 
transformation of "the time scale so that in the new time scale the data 
run is 2 seconds long, a process illustrated, in example II. 


Example II 


The first example served to illustrate the application of the method 
to an equation of-the form (l4), corresponding in the missile pitch- 
response problem to the case where only a(t) and 8(b) are measured. A 
problem involving equations, like ( 13 ), of the first order, corresponding 
to the case where q(t) is available in addition to a (t) and 6(t), will 
be illustrated now. 
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As In example I, the lift force will be assumed linear and the 
pitching moment nonlinear. The following parametric values were 
assumed: 


m 

= 2 

let = 1000 



V 

= 750 

M& = - 200 

► 

(29) 


= 100 

Mq = - 500 




The nonlinear M(a) is plotted as the solid curve in figure 3* It should 
be noted that in contrast to the first example , the pitching moment is 
unstable at a = 0 and highly nonlinear. 

The "test data" were manufactured by determining a pulse response of 
this missile on the REAC. Again, the control characteristics of the sys- 
tem will not be considered, so that only free oscillations are shown in 
figure 4. 

Merely to have a standard length of run, a 2- second interval was 
always selected for the calculation of the integrating factors r n (yv) 
(see the Appendix for the definition of these quantities). To illustrate 
the computation procedure for data runs of different lengths, a 1-Becond 
run will be considered in the present example. 


In order to use the integrating factors displayed in table III, it 
is necessary to make a preliminary transformation of the form 


t 



This transformation has the following effect on equations ( 13 ): 


■ ^ a - — a 2 - -i— a 3 - 2^ + q. = •— 5 

mV mV mV dT mV 

- a - “2 a? - !fea 3 - 2^^-^a + a£i = !^6 


(30) 


I y dr 


dT 


Recalling that for the problem under discussion, = 1^,, 112 = 13 = 0 , 
6 (t) = 0 , it can be seen from equations ( 30 ) that the equations to be 
solved by least squares for the parameters have the forms 
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- f a(T)y v (T)dT +2 J a(T ) dT + J q.(r)y v (r)6.r = 0 


4 


and 


- / a(T)y v (T)dT - ^ J a a (T)y v (T)dT - ^2 f a 3 (T)y v (r)dT + 

y 4 d y (d • L y <d 

2 1 r ■« ^ - - g/' *<^*<^ - 2 / 2 .(’> ^ « 


= 0 


The "data" of figure 4 are presented as functions of t in 
table III. The sums, which when multiplied by At = O.O 5 approximate 
the integrals in the last two equations, are given below table III. 

These sums have been listed in the appropriate places in table IV. With 
circled numbers referring to columns in table IV, it can be seen that 
the above two equations are equivalent to the following : 


- ^ © + 2 ® + © = 0 


mV 


© 

z y I y 


I 5 © + 2 

I y 


Mi 


-I® - 2 


Hence, 


la = mV 


2S©X® + 2©X© 

2 © 2 


or 


I<x = 1033 


= 0 


using the given values of m and V (eqs. (29) ), while the equations for 
the moment parameters are 



IZ ® 2 ‘l Z ®* ®*IZ 0*®-( 2 l)Z ®*® + IZ ®*®- 

- 2 £ © X ® 

%Z ®*®*IZ ®’ 4 IZ ®*©-( 2 l)Z ®*® 4 IZ ®*®- 

- e ^ © x © 


II 


^ ® x © +§ l 

\ © 2 -( s §)X © x ® + § £ ©x© = 




' 2 £ © X © > 

Mi \ 

‘ iy Z 

• © X ® -Y ) 

^ J y L 

? © X ® ^ 

-1 1 y L 

] © x © + (2 © 2 - Q £ © X ® = 




2 X ® X ® 

Mi V 

% z 

>*®*|I 

" © X © ] 

1 y L 

:®x«.(.s)z«x.^z®*- 




. 2 y © x @ 

u j 


When the sums displayed below table IV are inserted in the appropriate places, the resulting 
equations can be solved to yield 

UD 


NA.CA TN 3288 



2k 


MCA TN 3288 



521 

- 42 ^ 3 


2 ^ = 38.66 
I y 



26.7 


7 s = - 197,000 

J y 


Hence, using the value of* Iy given in equations (29), 


M z = 5-21X1Q 4 
Ms = - ^.23X10 3 


m 3 = - 1.97x10 7 


Hi = .1930 

Mq = - 2670 


(31) 


To begin our discussion of these values , consider L^. A compar- 
ison of the values given for by equations ( 29 ) and ( 3 l) shows that 
Ico has been found within 3*2 percent. The- nonlinear pitching moment 

M(cl) = (5.21Xl0 4 )a - ( ij-.23xlO s )cc 2 - (l.97Xl0 T )a 3 


has been superposed as the- dotted curve on the true moment curve in 
figure 3 . As can be seen, the agreement for this strongly nonlinear 
problem is excellent. Finally, the errors in the calculated values of 
the parameters 1%, and Mq are enormous, but- the reason for this Is well 
known and easily explained. Consider equations (13J which describe the 
motion. Eliminating q from these equations results in equation (l*0; 
since the lift ia linear , b x = b 2 = 0. The constant bo, on which Mq 
and Hi have their principal-effect, is easily interpreted physically 
as a measure of the damping in the system. The other gross aspects of— 
the response are relatively little affected by either Mq or.M<i. How- 
ever, the Important quantity in b Q is not Mq or Mi alone, but" is their 
sum, Mq + Mi* In other words, relatively large changes in Mq and Mi 
are possible without causing any great change in the motion, Just so long 
as their sum remains constant. Thus, one may not expec*tr-to find 
Mq and Mi accurately from an experiment such as this - only their sum 
may be relied upon. This is verified in the present example, since equa- 
tions ( 29 ) give the sum 


Mq + Hi = - 700 
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while equations (31) give 


Mq + M£ = - 740 


These two values differ by only 5.7 percent. 


It should be noted that the assumption, which has been made in 
this and the preceding example, of the linearity of the lift is not 
necessary. If it is suspected that the lift is linear but if no def- 
inite verification of this is available, a nonlinear form such as ( 12 ) 
can be assumed for the lift, and the coefficients L x , 1^, Lq, . . ., 
can be calculated. If the lift is indeed linear, it should turn out 
that 1^, Lq, . . ., are small. Some limited experience has shown that 
this method does work but that the errors in the calculated parameters 
are somewhat larger than when the correct form is* assumed. The reason 
for this is not known, but it appears that it may be associated with a 
tendency of the extraneous parameters Ig and L 3 to fit the lift curve 
to that corresponding to the original data, errors and all, at the 
expense of the smoothing operation which is necessary with this type of 
data and which is performed by the least squares process when the correct 
form is assumed. 


Example III 


We turn, finally, to a system described by a differential equation 
whose order is higher than the second. Since the higher order systems 
whose occurrence is most common appear to be those of the fourth order, 
we Bhall be concerned with such a system. In order to simplify the pres- 
entation, a linear system will be considered; there are no conceptual 
difficulties in the generalization to the nonlinear case. In addition. 

It will be assumed for simplicity that free oscillations are available 
for analysis. Thus, the system to be analyzed is assumed to be described 
by an equation of the form 


+ a 3 + a 2 + a x — + a Q x - 0 ( 32 ) 

dt 4 dt 3 dt 2 dt 

A solution of this equation was calculated over the interval from 0 
to 2 seconds, for the following values of the coefficients a^: 


a Q = 2544.9 
a 2 = 132.87 


a-L = 219.32 
a 3 = 2 . 000 
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The result, representing the free oscillations in response to some dis- 
turbance, is tabulated in table V and presented graphically as the solid 
curve in figure 5* The sums needed for the solution of the problem are 
displayed below table V and again in table VI. The least-squares equa- 
tions for the parameters are 


24.2035 a 0 - 16.7225 a x - 1098 . 19-32 + 1603.19 a Q 

- 16.7225 a Q + 855.538 a x + 664.495 ae - 82282.3 a 3 

- 1098.19 a Q + 664.495 + 77 B 98.7 ag - 94741.8 a 3 

1603.19 a Q - 82282.3 a x - 94741.8 ag + 852367 O a 3 


- 81036.4 
85777.1 
7533580 

- IO 5678 OO 


Solving these equations gives 


a Q = 3098.7 a! = 374^93 

aa = 141.29 ' ■ ■■ a 3 = 3*367 


It can be seen that these numbers are correct only to within orders o£- 
magnitude. On the other hand, it is not these coefficients which have 
direct, physical significance; rather, it - is the damping and frequency of 
each of the components making up the oscillation which are important. 

In order to find these numbers, the following equation was set up 


A 4 + 3-367 A 3 + lla .29 7\ s + 37^,93 A + 3098.7 = 0 
and solved to find the roots: 


- 0.046 ^10.7 i 

- 1.64 ± 4.96 i 


The true roots are obtained by solving the equation 


x 4 + 2.000 x 3 + 132.87 x 2 + 219.32 x + 2544^9 = 0 
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This gives 


X = 


± 10.5 i. 

- 1.00 ± 4.71 i 


Thus, the frequencies of the oscillation have been found quite accurately, 
as has the damping parameter of the undamped component. The only large 
effect of the errors in the coefficients a 0 , a-j_, a 2 , and a 3 is in the 
damping of one of the components. The apparent ill-conditioning of the 
problem with respect to this parameter is not too surprising, for after 
all, either the true or calculated value of this damping is large enough 
that the corresponding component of the motion is masked by the undamped 
component over a good part of the run. This may be seen best, perhaps, 
from figure 5 in which the solution of equation ( 32 ), using both the 
true and calculated values of the parameters, has been plotted. It can 
be seen that the two curves do not differ by very much, indicating that 
the fit could not be much improved. 


CONCLUDING REMARKS 


A general theory of the so-called "equations -of -motion" methods 
for the analysis of dynamical systems has been presented. It has been 
shown that, when looked at from a new point of view, all such methods 
can be generalized so as to apply to linear and nonlinear systems alike. 
Using this theory, it has also been shown how new methods can be devel- 
oped in order to satisfy the requirements of particular problems . 

One new method has been described in detail. In certain cases, it 
reduces to one which is very similar to the well-known. Fourier transform 
method (ref. 3) but in all cases has certain advantages over this latter 
method and over other methods heretofore used. Its superiority is based 
on two facts . First, there is the heavy dependence on the initial con- 
ditions which occurs when using most of the previously known equations - 
of-motion methods; this dependence is entirely eliminated in the new 
method. This superiority manifests itself particularly when systems of 
higher order than the second are considered. If, for definiteness, a 
fourth-order system is considered, before the Fourier transform method 
(for example) can be applied, it is necessary to evaluate the test data 
and their first three time derivatives at the initial point. Accurate 
evaluation of all these derivatives is practically impossible, however, 
with the type of data obtained from most aerodynamic experiments. 

The second fact upon which the superiority of the proposed method 
rests is that most of the equations -of-motion methods uBed to this time 
demand an infinitely long record for their rigorous application. For 
some years now, questions about the errors introduced into an analysis 
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of a system by the finite length of records available have been asked, 
but answers have not been offered. The second principal advantage of- 
the method described herein is that such questions are side-stepped com- 
pletely: There is no error at all from this source, since it-is assumed 

from the start that only a finite record is available. Because of this 
feature, the method avoids a further limitation of the- Fourier transform 
method (apparently the most accurate of all well.- known methods of this 
type), which cannot be applied at all to some systems (e.g., unstable 
ones) without time-consuming and sometimes ineffective special devices, 
since the Fourier .integrals of- the data simply do hot exist. 

The single exception to these remarks is the derivative method 
(refs. 2 and 3 ). The derivative method does noir^weight the initial con- 
ditions and does not depend on an infinite- interval for its application. 
In addition, the derivative method has in the past been considered as 
the only well-known method which applies to nonlinear as well as linear 
systems. (Other methods are described in references 2 and 10, but the 
derivative method appears' to be the only one with such general applica- 
bility as we are discussing here.) There are, however, a number of 
very serious objections to the derivative, method. First of all, there 
is the. inordinate amount of time and labor which must be expended in its 
application, principally because of the necessity for calculating time 
derivatives of the data. Second, and most important, is the question 
of accuracy. The accurate calculation of the derivatives needed for the 
method is most'difficult, and this calculation is a large source of 
error. Besides, even if the derivatives could be computed with the 
requisite accuracy, the derivative method appears often to lead to badly 
conditioned equations, as pointed out in reference 2 ; because of this, 
many problems have been found for which the derivative method has been 
shown to lead to extremely large errors. The method proposed herein is 
subject to none of these weaknesses. The time required for its applica- 
tion is far less than that needed for the derivative method; in addition, 
it appears to be well suited to machine computation. Naturally, deriv- 
atives need not be calculated, and the method shares the properties of 
the Fourier transform method which cause it to lead to fairly well- 
conditioned equations . 


Ames Aeronautical Laboratory 

National Advisory Committee for Aeronautics 
Moffett Field, Calif., July 1^- 195 ^ 



NACA TN 3288 


29 


APPENDIX 


NUMERICAL EVALUATION OF INTEGRALS OF THE FORM 



In 1928, Filon (ref. 8 - see also ref. 9 , pp. 67-72) published a 
generalization of Simpson* s rule for evaluation of integrals of the form 



x(t) sin cotdt 



x( t ) cos arfcdt 


where x(t) represents numerical data. Filon* s method, in contrast to 
Simpson*s, has the distinction of giving results whose errors are inde- 
pendent of the frequency co, depending only on how closely x(t) can be 
fitted to a sequence of parabolas. This method will be generalized to 
apply to integrals of the form 



x(t)y(t)dt 


where y(t) is known exactly (for application to the method described in 
the body of this report, y(t) is one of the method functions), while 
x(t) is given tabular ly. 


Suppose the interval (a,b) is divided into 2h equal parts by points 
t 0 = a < t x < . . . < t 2 k = b, where tn^ - t n = At = constant. Then, a 
formula of the form 



2h 

x(t)y(t)dt ffc At 

n=o 



r n (y)x(t n ) 


(33) 


will be sought, where the r n are constants which depend only on the 
function y. These constants will be determiried by the condition that 
the formula ( 33 ) will give the integral exactly in the cases where x(t) 
is a constant, a linear function or a quadratic function of t. Suppose 
first that the interval (a,b) is divided into two parts only by points 
t G , t x , and tg; since formula ( 33 ) is to be exact if x(t) = 1 , (t - t x ), 
or (t - t x ) 2 , we have 
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r G (y) + r x (y) + T a (y) = 2 y(t)at- 

to 

- (At)r 0 (y) + (At)r a (y) = ^ J' (t - t 1 )y(t)at- 


> (3^) 


(At) ! 


! r 0 (y) + (At) 2 r a (y) = J' (t - t 1 ) e y(t)atj 


Equations (34) can be solved for T 0f T 2 to obtain 

r ° (y) = / 2 (t ' tl)Sy(t)dt - f* (t ' tl)y(fr)dt 


r i(y) m EE f * y(t)dt - J^yS f a - (t - 


r a(y) - “ 7 Y 3 /" (t - t 1 ) £ y(t)dt + — i— /° (t - t^)y(t)dt 

2 (At ) s J to 2(At) 2 01 


>(35) 


t 0 




Now, if (a,b) is divided into 2h(>2) parts, the integral is written as 
follows : 

/ b r>t 2 nt4 

*(t)y(t)dt - J x(t)y(tr)dt + j x(t)y(t)dt + . . . + 


^t^h 


x(t)y(t)dt 

t2b-2 


(36) 


* 
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and equations (35) can be used to evaluate each, of the integrals on the 
right-hand side. Calling" 


Jp(y) - jj: / P+i y( t )^t 

■tp-i 

Kp(y) = — — '-g / (t - t p )y(t)at 

2 (At) J, 


tp-i 


1 nt p +1 

L p(y) = 3 / (t - t p ) 2 y(t)dt 

P ( At 'i l/. 


2 (At) 


bp-l 


we obtain from (33), (35 ^ and (36) that 


b 

^ x(t)y(t)dt « At ^ r n (y)x(t n ) 
a n=o 


where 


r 0 (y) = L x (y) - Kjy) 

r^.^y) = J 2p_i(y) - 2h= P -iiy)> p = i> 2, . . h 

r ap (y) = L ap _i(y) + K^>-i(y) + L ap+ 1 (y) - K apfi (y), 
P - 1, 2 , . . h - 1 

r 2h(y) = L sh..i(y) + K ah-i(y) 


> (37) 


It should he noted that if y(t) is identically unity. 


Jp(l) = 2 
Kp(l) = 0 
Lp(l) = 1/3 
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where the relations tp_x = tp - At, tp+x = tp + At have been used 
repeatedly. Hence, from equations (37) , 

r 0 (i) = 1/3 ; 

r ap-i(l) = ^/3j P - lj 2, . . . , h 

r 2 p (l) = 2 / 3 , p = 1 , 2 , . . h - 1 
r 2 h (i) = 1/3 - 

which exactly describes Simpson’s rule. 

Equations (37), with y(t) being chosen equal to the method func- 
tions, were used to calculate the numbers displayed, in tables I, III. 
and V. 
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TABLE V.- CALCULATION OF THE INTEGRALS NEEDED FOR EXAMPLE III - Concluded 


1 

72 

73 

74 

75 

76 

77 

78 

n 

—15 r n (y 13 ) 

r «(»u ( ‘>) 


■£“ r n(yze) 

& *W5»> 

^5 r n(yie) 

"1? r n(L (M ) 

0 

1.7865 

8.1426 

-0.0335 

-0.0693 

0.1514 

1.9516 

7.1188 

1 

.1117 

-19.8366 

.2580 

.8847 

1.7696 

-1.2097 

-19-9014 

2 

-7.931^ 

2.3892 

.5965 

-773^ 

-1.6842 

-6.7667 

11.5107 

3 

5.1605 

27.5582 

.9844 

-.9808 

-2.2435 

7.6917 

13.8720 

4 

4.0571 

-21.7377 

.0621 

-.5670 

1.5329 

-.0022 

-18.6294 

5 

-2.7607 

4.8180 

.0154 

0 

.9478 

0 

12.0588 

6 

3.5401 

• 5843 

.0621 

.5670 

1-5329 

.0022 

-18.6294 

7 

-7.3382 

-18.6211 

.9844 

.9808 

-2.2435 

-7.6917 

13.8720 
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-2.3484 

27.5345 

.5965 

-.7734 

-1.6842 

6.7667 

11.5107 
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9.3367 

-4.4593 

.2580 

-.8847 

1.7696 

1.2097 

-19.9014 

10 
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-12.4980 

-.0671 
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.3027 

0 

14.2374 

11 

1.1638 

13.3047 

.2580 

.8847 

1.7696 

-1.2097 

-19-9014 

12 

-3585 

-21.3858 

.5965 

• 7734 
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11.5107 

13 

-8.6181 

13.1777 
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-.9808 

-2.2435 

7.6917 

13.8720 

14 
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17.8633 
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1.5329 
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3-.65S2 - 
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0 

12.0588 

16 
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.5670 
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17 
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19 
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20 
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21 

-2.4995 

-8.2646 

.2580 

.8847 

1.7696 

-1.2097 

-19.9014 

22 

2.4454 

11.8676 

.5965 

.7734 

-1.6842 

-6.7667 

11.5107 

23 

-3.6582 

-24.2143 

.9844 

-.9808 

-2.2435 

7.6917 

13.8720 

2li 

-6.2166 

17.8633 

.0621 

-5670 

1.5329 

-.0022 

-18.6294 

25 

8.6181 

13.1777 

.0154 

0 

.9478 

0 

12.0588 

26 

.3585 

-21.3858 

.0621 

.5670 

1.5329 

.0022 

-18.6294 

27 

-1.1638 

13.3047 

.9844 

.9808 

-2.2435 

-7.6917 

13.8720 

28 

2.5377 

-12.4980 

• 5965 

-.7734 

-1.6842 

6.7667 

11.5107 

29 

-9.3367 

-4.4593 

.2580 

-.8847 

1.769 6 

1.2097 

-19.9014 

30 

2.3484 

27.5345 

—0671 

0 

.3027 

0 

14.2374 

31 

7.3382 

-18.6211 

.2580 

.8847 

1.7696 

-1.2097 

-19-9014 

32 

-3.5401 

.5843 

.5965 

.7734 

-1.6842 

-6.7667 

11.5107 

33 

2.7607 

4.8180 

.9844 

-.9808 

-2.2435 

7.6917 

13.8720 

3k 

-4.0571 

-21.7377 

.0621 

-5670 

1.5329 

-.0023 

-18.6294 

35 

-5.1605 
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0 
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0 
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36 
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1.5329 

.0022 

-18.6294 

37 

-.1117 

-19.8366 
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.9808 

-2.2435 

-7.6917 

L3.872O 

38 

-1.7865 

8.1426 

.5965 

-*7734 

-1.6842 

6.7667 

11.5107 

39 

0 

0 

.2580 

-.8847 

1.7696 

1.2097 

-19.9014 

40 

0 

0 

-.0335 

.0693 

.1514 

-1.9516 

7.1188 

O 

= -.055 

.169 . 

-.709 

.2724 

-.0128 

-.114 

.235 





333*933- 

093*^9- 

£130*3- 

3l6£’T- 


C36*39£- 

9*ii£*3- 

6t£1*£ 

900*T£ 

0T3TT3 

£rri* 

616* l£i~ 

£8£3*£- 

T3*l3*1 

8££*l£ 

9889' IS 

8W£‘ 

£9T*3^- 

9£63*3I- 

03H*£ 

oli*6£ 

6931*33 

9193* 

119*3911“ 

9958’ £<r 

£969*01 

osrsgn- 

6£l0'£33 

£g£*?*Tt 

16r£l6l 

H£t*t£t 

0663*91- 

90f9U 

393£*3t- 

1668*9“ 

s£9‘«i- 

89l£*19 

6££o*oi- 

llr£s 

0£t£*££- 

3£lO*l- 

£01*31- 

lgo£*i3- 

9£30*I- 

< 3 > x © 

© X © 

© x © 

3t 

Tt 

01 


££s* 

69T* 

66i* 

6l0'- 

8l£* 

w- 

906* 

£6r- 

lM*1 

100*13- 
033*13“ 
100*9T 
61£*6£- 
£69*001 
903* 601 


OT^OT- 

9*699£3ee 

5*9 /tfEKZ- 
a£l*TlH6 
931*86811 
Torui£g 
3£E* 38330- 


@ x @3 
@ X @2 

s@3 
@ x @ 3 

(§T) X (oft 2 


E9n6rw 

TX&C'QZQ 
0££*9£ OT0 
CE6l*£09X- 
3991 '8601- 
10£33l‘9r 
6o££03*13 



@ X © 2 

@ x © 2 

@ X © 3 
3© 3 


ITT*- 

03TO*- 

13l3* 

E6ol*- 

313^*9^6*13 

1io*t**i96i 

66££*o££l 

l£l6*l£l 

*l99£*3T 

££o*- 

6ttq*- 

0I3T 

6 £ 69 - 

6£6£’9l6lT 

3*i96*££t 

El 2 £*TC 

£13*- 

8030*- 

l6££* 

3901*- 

lT£*i*lI9nl 

T*t 6 £* 6 s£t 

l306*Q3t 

9 £ 66 *OT 

lao* 

9361* 

£6£o* 

9^8- 

2£T8’£303I 

909£-8*tTI 

£399*601 

03l*r*ot 

9i£*- 

36£o*- 

io£l* 

£ 969 - 

£ 9 £l‘ 068 l 

£69I*1£8 

1938*88 

913*f*6 

a£o* 

8193* 

ffnio* 

9106 *- 

8£01 * 661£ 

I3l£**i99 

*rt£l*9l 
0^89 *19 

9931*8 

xgi’T- 

E 66 t— 

6369* 

H99*- 

93»i0*£08£ 

T£iir*t 8 n 

0 *l£ 8 *l 

9 TC* 

£991* 

£l£o* 

L9£6*~ 

* 169*1 *£l£3 

609 s* 0 *|£ 

90£1*91 

£ 196*9 

2.89*1- 

60TT*T- 

130l*T 

£3 &r- 

*?£*(£ * 8 ££i 

S0£o*9*[3 

10l1*6E 

3£83’9 

£ 16 * 9 - 

oliE’i 

8990*3 

30£8*3- 

6089 '^6 

986£*o1t 

8 ££ 1 *o£ 

0 *rl£*£ 

£ 18 * 8 * 

a£o6*£ 

3£98*£“ 

192**3- 

6£Er£6r 

3919* lot 

9903*33 

ieil *1 

li 9 *01 

£661*- 

Ol£l*I- 

£ 8 £ 8 ‘- 

a£i8*l£3 

I6££’09 

£i3*i*£t 

0l36*£ 

000*13 

16t£*9 

6003*£- 

9T£i*3- 

I 601* 16 

£ 900 *T£ 

9698*6 

’ 9W*£ 

319** 

06 £ 1 *£~ 

£ll6’s- 

398T*- 

9£ot**t£ 

0 £TT*tT 

ooie*£ 

9911*3 

196*01- 

9^90* ii“ 

8l£9’0- 

13*10’ I 

1890*9 

8 £l 0 *£ 

12-91* s 

901S*I 

« (0) «r Jr w 3 * 
2 

r V7— 
T 

q-P^fx r q.7 

3 J 1 


A 



Aw 




saavnbs iswi is in aiawvxa so saaiawvavd aoai jd mirnmm ms aaaaaN swns -*ia srcera 


99s£ Hi vdvh 


































ft-lb 


NACA TN 3288 


^7 



Figure 1.* The nonlinear pitching ^moment curve for example I. 
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Figure k,- 


"Test data" for example II 





Figure 5.- "Test data" for example III 




